#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use cluster_from_pairs;
use parse_bl2seq;
use chromosomes;
use DBI;
use DBD::Pg;
use linewrap;
use Spreadsheet::WriteExcel::Big;


(my $my_dir = $FindBin::Bin) =~ s/(.*)\/$/$1/;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    xls_protocol_synteny_histograms

	Description:

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------


	use constant PAIR_NAME		=> 0;
	use constant PAIR_NAMES		=> 1;
	use constant PAIR_IDS		=> 2;








#_________________________________________________________________________________________

#	get_xls_formats

#_________________________________________________________________________________________
sub get_xls_formats(\%$)
{
	my ($formats,  $workbook) = @_;

	# bold
	$formats->{bold} = $workbook->addformat();
	$formats->{bold}->set_bold();
	
	# blue and bold
	$formats->{blue_bold} = $workbook->addformat();
	$formats->{blue_bold}->set_bold();
	$formats->{blue_bold}->set_color('blue');
	
	# red and bold on yellow background
	$formats->{red_bold_yellow} = $workbook->addformat();
	$formats->{red_bold_yellow}->set_bold();
	$formats->{red_bold_yellow}->set_color('red');
	$formats->{red_bold_yellow}->set_bg_color('yellow');
	
	# three significant digits
	$formats->{double_digit} = $workbook->addformat();
	$formats->{double_digit}->set_num_format("0.00");
	
	# percentage
	$formats->{percent} = $workbook->addformat();
	$formats->{percent}->set_num_format("0%");
	
	# percentage
	$formats->{blue_bold_percent} = $workbook->addformat();
	$formats->{blue_bold_percent}->set_bold();
	$formats->{blue_bold_percent}->set_color('blue');
	$formats->{blue_bold_percent}->set_num_format("0%");

	# normal
	$formats->{black} = $workbook->addformat();
}

#_________________________________________________________________________________________

#	get_histogram

	# calculate counts, quartiles and percentile histograms
	# output format is counts in one array, 
	#			quartiles and histograms in another array but separated by 2 ''

#_________________________________________________________________________________________
sub get_histogram(@)
{
	my @data = sort {$a <=> $b} grep {defined $_} @_;
	my $count_undefined = scalar grep {!defined $_} @_;
	my $count = scalar @data;
	my @indices = (@data * 0.25, @data * 0.5, @data * 0.75);
	my @quartiles = @data[@indices];
	
	@indices = (0..100);
	$_ = int($_ * @data / 100) for @indices;
	$indices[-1] = @data - 1;
	$indices[0] = 0;
	my @percentiles = @data[@indices];

	return [[$count, $count_undefined], [@quartiles, '', '', @percentiles]];
}



#_________________________________________________________________________________________

#	do_select_subset

#		get count or histogram of data divided by all, protocol1 only , protocol2 only, 
#				common or all protocol1, all protocol2

#_________________________________________________________________________________________
sub do_select_subset_by_protocol_id($$$$$$)
{
	# new
	my ($dbh, $return_data, $protocol_id_name, $protocol_id, $what, 
		$table) = @_;

	my $sql_cmd = <<"PL/SQLCMD";
        SELECT 
            $what, 
            oso.species
        FROM
            $table
            NATURAL JOIN orthologs.real_genes
            JOIN orthologs.one_or_many_orthologs oso  USING (gene_id) 
        WHERE
            oso.$protocol_id_name & $protocol_id = $protocol_id AND
            a.protocol_id = $protocol_id;
PL/SQLCMD
	#print STDERR "$sql_cmd\n";
	my $data = $dbh->selectall_arrayref($sql_cmd);


	#orthologues unique to this protocol from both species
	{	
		my @test_data = map{$_->[0]} @$data;
		push(@$return_data, get_histogram(@test_data));						
	}

	#all orthologues from species 1 unique to this protocol
	{	
		my @test_data = map{$_->[0]} grep {$_->[1] eq $ortho_name1} @$data;
		push(@$return_data, get_histogram(@test_data));						
	}

	#all orthologues from species 2 unique to this protocol
	{	
		my @test_data = map{$_->[0]} grep {$_->[1] eq $ortho_name2} @$data;
		push(@$return_data, get_histogram(@test_data));						
	}


}
sub do_select_subset($$$$$$$$)
{
	my ($dbh, $return_data, $names, $protocol_name, 
		$what, $table, $protocol_id_name, $id) = @_;

	push (@$names, 	"" ,
					"$protocol_name",
					"$protocol_name unique $ortho_name1",
					"$protocol_name unique $ortho_name2",
					"Percentage");


	do_select_subset_by_protocol_id($dbh, $return_data, $protocol_id_name, $id, $what,
									$table);
	
	# null column at end
	push(@$return_data, [[''], ['']], [[''], ['']]);						

}

#_________________________________________________________________________________________

#	get_histograms_for_all_relationship_types

#		get count or histogram of data for orthologs, 1:1 orthologues, manys, paralogues

#_________________________________________________________________________________________
sub get_histograms_for_all_relationship_types($\@\@\@$$$$)
{
	my ($dbh, 
		$data, $names, # return data and names of columns
		$block_names,
		$protocol_name,
        $what, 
		$table,
		$protocol_id) = @_;

	# select orthologues
	print STDERR "\tGet orthologues data\n";
	do_select_subset($dbh, 
					$data, 
					$names, 
					$protocol_name, 
					$what, 
					$table,
					'protocol_id', 
					$protocol_id);
	push(@$block_names, "orthologues");

	# select 1:1 orthologues
	print STDERR "\tGet 1:1 orthologues data\n";
	do_select_subset($dbh, 
					$data, 
					$names, 
					$protocol_name, 
					$what, 
					$table,
					 'protocol_id_1s', 
					$protocol_id);
	push(@$block_names, "1:1 orthologues");

	# select manys orthologues
	print STDERR "\tGet manys orthologues data\n";
	do_select_subset($dbh, 
					$data, 
					$names, 
					$protocol_name, 
					$what, 
					$table,
					 'protocol_id_manys', 
					$protocol_id);
	push(@$block_names, "manys orthologues");
	return 3;
}

#_________________________________________________________________________________________

#	write_range_data_to_xls

#		send histogram etc to xls

#_________________________________________________________________________________________
sub write_range_data_to_xls($$$\@\@\@$)
{
	my ($worksheet, $format, $cnt_blocks, $block_names, $row_names, $data, $data_format) = @_;

	# write data
	{
		my $col = 1;
		for $data(@$data)
		{
			# counts
			$worksheet -> write_col(2, $col, $data->[0]);
			#quartiles and percentiles
			$worksheet -> write_col(4, $col, $data->[1], $data_format);
			++$col;
		}
	}

	# write percentages (Y column) every 14
	my @range = map {$_ / 100} (0..100);
	$worksheet->write_col(9, $_ * 5 - 1, \@range, $format->{blue_bold_percent}) 
																for (1.. $cnt_blocks);

	# write block names
	$worksheet->write(0, $_ * 5 + 1, $block_names->[$_], $format->{red_bold_yellow}) 
															for (0.. ($cnt_blocks - 1));
	# write header row and column
	{
		my $col = 0;
		for (@$row_names)
		{
			if ($_)
			{
				$worksheet->write(1, $col, $_, $format->{blue_bold});
				$worksheet->write(8, $col, $_, $format->{blue_bold});
			}
			++$col;
		}
	}
	$worksheet -> write_col(2, 0, 
							['count', 'count missing', '="25%"', "median", '="75%"'], 
							$format->{blue_bold});
}

#_________________________________________________________________________________________

#	get_blocksize_histogram

#		get histogram of block sizes

#_________________________________________________________________________________________
sub get_blocksize_histogram($$\%$)
{
	print STDERR "Get block size histogram...\n";
	my ($dbh, $workbook, $format, $protocol) = @_;
	my ($name_prefix, $id, $name) = @$protocol;

	# retrieved data
	my (@data, @names, @block_names);

	my $cnt_blocks = get_histograms_for_all_relationship_types
								($dbh, @data, @names, @block_names,
								 $name,
								 'block_size', 
								 "orthologs.micro_syntenic_blocks a JOIN\n\t\t".
								 'orthologs.micro_syntenic_block_sizes '.
								 'USING (micro_syntenic_block_id, species, protocol_id)',
								 $id);
	my $worksheet = $workbook->addworksheet($name_prefix . " block size");
	write_range_data_to_xls($worksheet, $format, $cnt_blocks, 
						@block_names, @names, @data, 
						undef);
}


#_________________________________________________________________________________________

#	get_closest_neighbour_histogram

#		get histogram of closest neighbour

#_________________________________________________________________________________________
sub get_closest_neighbour_histogram($$\%$)
{
	print STDERR "Get closest neighbour distance histogram...\n";
	my ($dbh, $workbook, $format, $protocol) = @_;
	my ($name_prefix, $id, $name) = @$protocol;

	# retrieved data
	my (@data, @names, @block_names);

	my $cnt_blocks = get_histograms_for_all_relationship_types
								($dbh, @data, @names, @block_names,
								 $name,
								 'closest_neighbour_dist',
								 'orthologs.syntenic_distances a ',
								 $id);
	my $worksheet = $workbook->addworksheet($name_prefix . " closest neighbour");
	write_range_data_to_xls($worksheet, $format, $cnt_blocks, 
						@block_names, @names, @data, 
						undef);
}




#_________________________________________________________________________________________

#	get_triplet_histogram

#		get histogram of triplets

#_________________________________________________________________________________________
sub get_triplet_histogram($$\%$)
{
	print STDERR "Get triplets histogram...\n";
	my ($dbh, $workbook, $format, $protocol) = @_;
	my ($name_prefix, $id, $name) = @$protocol;

	# retrieved data
	my (@data, @names, @block_names);

	my $cnt_blocks = get_histograms_for_all_relationship_types
								($dbh, @data, @names, @block_names,
								 $name,
								 'other_dir_neighbour_dist', 
								 'orthologs.syntenic_distances a ',
								 $id);
	my $worksheet = $workbook->addworksheet($name_prefix . " triplet");
	write_range_data_to_xls($worksheet, $format, $cnt_blocks, 
						@block_names, @names, @data, 
						undef);
}

#_________________________________________________________________________________________

#	get_syntenic_distance_counts

#_________________________________________________________________________________________
sub get_syntenic_distance_counts($$$$$$$)
{
	my ($dbh, $total_id, $protocol_id, $protocol_id_name, $unique, $where, $table) = @_;
	my $where_unique = $unique ? "AND oso.protocol_id & $total_id = $protocol_id" : "";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
		COUNT(*),
		oso.species
	FROM
		$table
		NATURAL JOIN orthologs.real_genes
		JOIN orthologs.one_or_many_orthologs oso  USING (gene_id) 
	WHERE
		oso.$protocol_id_name & $protocol_id <> 0 AND
		a.protocol_id = $protocol_id
		$where_unique
		$where
	GROUP BY 
		oso.species;
PL/SQLCMD
#	print STDERR "$sql_cmd\n";
#	return (100,100,100);
	my $data = $dbh->selectall_arrayref($sql_cmd);
	my ($total_cnt, $cnt1,  $cnt2) = (0,0,0);
	for (@$data)
	{
		if ($_->[1] eq $ortho_name1)
		{
			$cnt1 +=$_->[0];
		}
		else
		{
			$cnt2 +=$_->[0];
		}
		$total_cnt +=$_->[0];
	}
	return ($total_cnt, $cnt1,  $cnt2);
}


sub get_closest_neighbour_summary($$\%\@)
{
	print STDERR "get_closest_neighbour_summary...\n";
	my ($dbh, $workbook, $format, $protocols) = @_;

	my $worksheet = $workbook->addworksheet("Closest Neighbour Summary");
	
	my $row_ortholog = 1;
	my @ortholog_types = (["orthologs", "protocol_id"],
						 ["1:1 orthologs", "protocol_id_1s"],
						 ["manys orthologs", "protocol_id_manys"]);
	
	for my $ortholog_type(@ortholog_types)
	{
		my ($ortholog_name, $protocol_id_name) = @$ortholog_type;

		# ortholog_type
		$worksheet->write($row_ortholog, 1, $ortholog_name, $format->{blue_bold});


		
		# column header
		my @col_header = ("all genes", "< 5", "< 50", "same chromosome", "other chromosome");
		my $offset1 = 8;
		my $offset2 = 16;
		$worksheet->write_col($row_ortholog + 3, 1, \@col_header, $format->{red_bold_yellow});
		$worksheet->write_col($row_ortholog + 3 + $offset1, 1, \@col_header, $format->{red_bold_yellow});
		$worksheet->write_col($row_ortholog + 3 + $offset2, 1, \@col_header, $format->{red_bold_yellow});
		#species name
		$worksheet->write($row_ortholog + 1 + $offset1, 1, $ortho_name1, $format->{red_bold_yellow});
		$worksheet->write($row_ortholog + 1 + $offset2, 1, $ortho_name2, $format->{red_bold_yellow});
		
	
		#
		# get total id
		#
		my $total_id = 0;
		for my $protocol(@$protocols)
		{
			$total_id += $protocol->[1];
		}
	
		#
		# row headers
		#
		my $col = 2;
		my $unique_col_offset = @$protocols * 2;
		for my $protocol(@$protocols)
		{
			my ($protocol_name) = $protocol->[0];
			for (0..2)
			{
				$worksheet->write($row_ortholog + 2 + $_ * $offset1, $col, $protocol_name, $format->{blue_bold});
				$worksheet->write($row_ortholog + 2 + $_ * $offset1, $col + $unique_col_offset, $protocol_name, $format->{blue_bold});
			}
			$col += 2;
		}
		$worksheet->write($row_ortholog + 1, 2, "All genes", $format->{blue_bold});
		$worksheet->write($row_ortholog + 1, 2 + $unique_col_offset, "Unique genes", $format->{blue_bold});



		$col = 0;
		{
			#
			# for each protocol
			#
			for my $protocol(@$protocols)
			{
				$col += 2;
				my ($protocol_name, $protocol_id) = @$protocol;
		
				# all counts unique or otherwise
				my (@counts_all) = 
					get_syntenic_distance_counts($dbh, $total_id, $protocol_id, 
												 $protocol_id_name, 0, "", 
												 " orthologs.syntenic_distances a ");
				my (@counts_all_unique) = 
					get_syntenic_distance_counts($dbh, $total_id, $protocol_id, 
												 $protocol_id_name, 1, "", 
												 " orthologs.syntenic_distances a ");
	
				# write counts by both and individual species
				my $row = 3;
				for (0 .. 2)
				{
					$worksheet->write($row_ortholog + $row + $_ * $offset1, $col,	
									  $counts_all[$_]);
					$worksheet->write($row_ortholog + $row + $_ * $offset1, $col + 1, 	
									  "1.0", $format->{percent});
					$worksheet->write($row_ortholog + $row + $_ * $offset1, $col + $unique_col_offset, 
									  $counts_all_unique[$_]);
					$worksheet->write($row_ortholog + $row + $_ * $offset1, $col + $unique_col_offset + 1, 
									  "1.0", $format->{percent});
				}
	
				# do same for different conditions
				for my $condition('< 5', '< 50', '< 20000', ' = 20000')
				{
					++$row;
					my $where = ' AND closest_neighbour_dist '. $condition . ' ';
					# other_dir_neighbour_dist
	
					my (@counts) = 
						get_syntenic_distance_counts($dbh, $total_id, $protocol_id, 
													 $protocol_id_name, 0, $where, 
													 " orthologs.syntenic_distances a ");
					my (@counts_unique) = 
						get_syntenic_distance_counts($dbh, $total_id, $protocol_id, 
													 $protocol_id_name, 1, $where, 
													 " orthologs.syntenic_distances a ");
					for (0 .. 2)
					{
						$worksheet->write($row_ortholog + $row + $_ * $offset1, $col,	
										  $counts[$_]);
						$worksheet->write($row_ortholog + $row + $_ * $offset1, $col + 1, 	
										  $counts[$_] / $counts_all[$_], 
												$format->{percent});
						$worksheet->write($row_ortholog + $row + $_ * $offset1, $col + $unique_col_offset, 
										  $counts_unique[$_]);
						$worksheet->write($row_ortholog + $row + $_ * $offset1, $col + $unique_col_offset + 1, 	
										  $counts_unique[$_] / $counts_all_unique[$_], 
												$format->{percent});
					}
				}
				
			}
		}
		$row_ortholog += 27;
	}
}




#_________________________________________________________________________________________

#	get_triplet_histogram

#		get histogram of triplets

#_________________________________________________________________________________________
sub get_macro_histogram($$\%\@)
{
	print STDERR "Get macro histogram...\n";
	my ($dbh, $workbook, $format, $protocols) = @_;

	my $col = 1;
	my $worksheet = $workbook->addworksheet("Macro block histogram");
	for my $species($ortho_name1, $ortho_name2)
	{
		# species on first row
		$worksheet->write(1, $col, $species. ' macro sizes', $format->{red_bold_yellow});
		for my $protocol(@$protocols)
		{
			my ($protocol_name, $protocol_id) = @$protocol;
			my $sql_cmd = <<"PL/SQLCMD";
			SELECT 
				block_size, 
				cummulative_gene_fraction
			FROM
				orthologs.macro_syntenic_block_size_histograms
			WHERE
				protocol_id = $protocol_id AND
				species = '$species';
PL/SQLCMD
			my $data = $dbh->selectall_arrayref($sql_cmd);	

			# column name on second row
			$worksheet->write(2, $col, $protocol_name, $format->{blue_bold});
			$worksheet->write(2, $col +1, "% of Genes", $format->{blue_bold});

			my @block_sizes = map {$_->[0]} @$data;
			my @fraction 	= map {$_->[1]} @$data;
			$worksheet->write_col(3, $col, \@block_sizes);
			$worksheet->write_col(3, $col + 1, \@fraction, $format->{percent});
			$col +=2;
		}
		$col +=1;
	}
}




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Get data from panda
#
#
my $dbh = connect_to_panda();

my $workbook  = Spreadsheet::WriteExcel::Big->new("$FindBin::Bin/../paper/block_sizes_syntenic_distance.xls");
die "Problems creating new Excel file: $!" unless defined $workbook;

my %formats;
get_xls_formats(%formats, $workbook);

$dbh->do ("set work_mem=160384");


	#
	# get protocol_ids and comparison pairs
	#
	my %protocol_names_to_ids;
	$protocol_names_to_ids{$_->[0]} = $_->[1] for (get_orthologs_protocols());
	die unless exists $protocol_names_to_ids{'phyop'};

	my @protocols = (
					["phyop", $protocol_names_to_ids{'phyop'}, "phyop"]
				);

	for my $protocol(@protocols)
	{
		print "\n\n$protocol->[0]...\n";
		get_blocksize_histogram($dbh, $workbook, %formats, $protocol);
		get_closest_neighbour_histogram($dbh, $workbook, %formats, $protocol);
		get_triplet_histogram($dbh, $workbook, %formats, $protocol);
	}

	my @selected_protocols = (
						   ["phyop", $protocol_names_to_ids{'phyop'}],
							);
	get_macro_histogram($dbh, $workbook, %formats, @selected_protocols);

	shift(@selected_protocols);
	get_closest_neighbour_summary($dbh, $workbook, %formats, @selected_protocols);
	#get_triplet_summary($dbh, $workbook, %formats, @selected_protocols);


	

$workbook->close();
print STDERR "Finished.\n";


